#'###########################################
#'                                          #
#'  ------------------------------------    #
#'  |     COIN     __                  |    #
#'  |            >(' )                 |    #
#'  |              )/                  |    #
#'  |             /  `----/            |    #
#'  |             \  ~=- /             |    #
#'  |   ~^~^~^~^~^~^~^~^~^~^~^~^~^~^   |    #
#'  ------------------------------------    #
#'                                          #
#'###########################################

setClassUnion("CharOrNum",   c("character", "integer"))

HRCompareParams <- setClass("HRCompareParams", representation(
				iterations="integer",
				mcpPercent="integer",
				kernelPercent="integer",
				kernelh="CharOrNum",
				kernelCellsize="integer",
				kernelExtent="numeric",
        kernelExtentCell="integer",
				doKernel="logical",
				doMcp="logical",
				cases="list"))



#'#
#'# Compare Argos and GPS HR with iterations on the generation of the GPS HRs for
#'# mobile data
#'# 
compareHRMobile <- function(data, iterations = 100L,
		                        resultsPath=".", writeResults=FALSE, name="", 
		                        mcpPercent=100L, kernelPercent=95L, kernelh="href",
                            kernelCellsize=250L, kernelExtent=0.3,
                            kernelExtentCell=10L) {
  
  	print("Comparing Home ranges")
  	
  	hrCompareParams <- HRCompareParams(iterations, mcpPercent, kernelPercent, 
  			kernelh, kernelCellsize, kernelExtent, kernelExtentCell, "Mobile")
  	
  	#' Comparison of HR by deployment
    tmp <- lapply(data, compareHomeRanges, hrCompareParams)

    res <- rbindlist(tmp)
    
  	if (writeResults){
  		fileName <- paste("HR_compare_", name, ".csv", sep="")
  		write.csv(res, file.path(resultsPath, fileName))
  	}
  	return (res)
}



#'#
#'# Generates Argos HR and compare it with randomly generated GPS HR
#'# 
compareHomeRanges <- function (data, hrCompareParams) {

    
    print(paste("Computing for deployment", data$deployment))
    
		#' if we have less than 5 locs, we can't calculate a HR
		if (is.null(data$compared)) {
			#'# TODO A refaire
			#res <- as.list(as.double(c(n, rep(NA, length(HRCOMPARE_COLNAME) - 1))))
			#names(res) <- HRCOMPARE_COLNAME
			#return (res)
      return(NULL)
		}
    
    
    n <- nrow(data$compared@coords)

		#' Transform into SpatialPoints
		#' comparedLocs <- SpatialPoints(data$compared)
		
		#' Argos HR via MCP. No need to do it more than once
    comparedMcp <- NULL
		if (hrCompareParams@doMcp) {
			comparedMcp <- lapply(hrCompareParams@mcpPercent,
                         function(percent, locs){
                           mcp(locs, percent=percent, unout="m2")
                         }, locs=data$compared)
			names(comparedMcp) <- hrCompareParams@mcpPercent
		} 
    
		#' Argos HR via Kernel
		kernelGrid <- NULL
		argosKernel <- NULL
		if (hrCompareParams@doKernel) {
			kernelGrid <- makeKernelGrid(c(data$ref@coords[, 1], data$compared@coords[, 1]),
                                   c(data$ref@coords[, 2], data$compared@coords[, 2]),
                                   hrCompareParams@kernelCellsize,
                                   hrCompareParams@kernelExtent)
			comparedKernel <- kernelUD(data$compared,
                              h=hrCompareParams@kernelh,
                              grid=kernelGrid)
		}
		
		#' Multithreading si possible pour faire les iterations
		 if (exists("ARGOS_CLUSTER")) {
		    tmp <- parLapply(ARGOS_CLUSTER, 1:hrCompareParams@iterations,
                         compareWithReference, data=data, n=n,
                         comparedMcp=comparedMcp, comparedKernel=comparedKernel,
                         kernelGrid=kernelGrid, hrCompareParams=hrCompareParams)
		} else {
			  tmp <- lapply(1:hrCompareParams@iterations, compareWithReference, data=data,
                      n=n, comparedMcp=comparedMcp, comparedKernel=comparedKernel,
                      kernelGrid=kernelGrid, hrCompareParams=hrCompareParams)
	  }
		
    res <- do.call(rbind, tmp)
    res$deployment <- data$deployment

    return(res)
}



#'#
#'# Performs an iteration of GPSS HR generation and comparison
#'# 
compareWithReference <- function(iteration, data, n, comparedMcp, comparedKernel,
                           kernelGrid, hrCompareParams) {
	
		print(paste("Iteration:", iteration))
	
		#' Sample a number of gps points equals to the number of Argos locs
		rndIndexes <- sample(nrow(data$ref@coords), n)
    data$ref <- data$ref[rndIndexes, ]

    #' Kernel
    refKernel <- NULL
		if (hrCompareParams@doKernel) {
        refKernel <- kernelUD(data$ref, h=hrCompareParams@kernelh, grid=kernelGrid)
		}
		
		tmp <- lapply(hrCompareParams@cases, getHRStatistics, data=data, n=n,
                  comparedMcp=comparedMcp, comparedKernel=comparedKernel,
                  refKernel=refKernel)
    
    res <- do.call(rbind, tmp)
		
    return(res)
}



#'#
#'# Creates an object with parameter for HR comparisons
#'# 
HRCompareParams <- function (iterations = 100L, mcpPercent=100L, kernelPercent=99L, 
                             kernelh="href", kernelCellsize=250L, kernelExtent=0.3,
                             kernelExtentCell=10L, type) {
    hrParams <- new("HRCompareParams")
    hrParams@iterations <- iterations
    hrParams@doMcp <- (length(mcpPercent) > 0)
    hrParams@doKernel <- (length(kernelPercent) > 0)
    if (hrParams@doMcp) {
      hrParams@mcpPercent <- mcpPercent
    }
    if (hrParams@doKernel) {
      hrParams@kernelCellsize <- kernelCellsize
      hrParams@kernelPercent <- kernelPercent
      hrParams@kernelh <- kernelh
      hrParams@kernelExtent <- kernelExtent
      hrParams@kernelExtentCell <- kernelExtentCell
    }
    hrParams@cases <- createCases(mcpPercent, kernelPercent, type)
    return (hrParams)
}


#'#
#'# Returns a list of home ranges to compare based on method and percentage
#'# 
createCases <- function(mcpPercent, kernelPercent, type) {
    mcp <- as.list(mcpPercent)
    for (i in seq(along=mcp)) {
      class(mcp[[i]]) <- paste("mcpPercent", type, sep="")
    }
    kernel <- as.list(kernelPercent)
    for (j in seq(along=kernel)) {
      class(kernel[[j]]) <- paste("kernelPercent", type, sep="")
    }
    return (c(mcp, kernel))
}


#'#
#'# Generic method for home range statistics
#'# 
getHRStatistics <- function(...){
	  UseMethod("getHRStatistics")
}


#'#
#'# Compute statistics for the MCP home ranges
#'# 
getHRStatistics.mcpPercentMobile <- function(percent, n, data, comparedMcp, ...) {

    #' Home Ranges
    refHR <- mcp(data$ref, percent=percent, unout="m2")
    comparedHR <- comparedMcp[[as.character(percent)]]
   
    
    #' Areas
    refArea      <- refHR@data$area
    comparedArea <- comparedHR@data$area
    
    ## print(data$ref@coords)
     ## plot(data$compared, xlim=c(min(data$compared@coords[,"x"], data$ref@coords[,"x"]), max(data$compared@coords[,"x"], data$ref@coords[,"x"])), 
     ##     ylim=c(min(data$compared@coords[,"y"], data$ref@coords[,"y"]), max(data$compared@coords[,"y"], data$ref@coords[,"y"])))
     ## points(data$ref, col="green")
     ## plot(comparedHR, add=T)
     ##  plot(refHR, border="blue", add=T)

    
    class(percent) <- NULL
    return (getHRStats(percent, "mcp", n, refHR, refArea, data$compared@coords,
            comparedHR, comparedArea))
}


#'#
#'# Compute statistics for the kernel home ranges
#'# 
getHRStatistics.kernelPercentMobile <- function(percent, data, n, refKernel, comparedKernel, ...) {

    #' Home Ranges
    refHR      <- getverticeshr(refKernel, percent=percent, unout="m2")
    comparedHR <- getverticeshr(comparedKernel, percent=percent, unout="m2")
    
    #' Areas
    refArea      <- refHR@data$area
    comparedArea <- comparedHR@data$area
    
    ##print(comparedKernel@h$h)
       ## image(comparedKernel)
       ## plot(comparedHR, add=T)
       ## plot(refHR, border="blue", add=T)
       ## points(data$compared)
       ## points(data$ref, col="green")
      
    class(percent) <- NULL
    
    return (getHRStats(percent, "kernel", n, refHR, refArea, data$compared@coords,
            comparedHR, comparedArea))
}



#'#
#'# Compute statistics for the home ranges
#'# 
getHRStats <- function(percent, method, n, refHR, refArea, comparedLocs, comparedHR,
    comparedArea) {

    #' Size ration Argos/GPS
    sizeRatio <- comparedArea / refArea

    #' Overlap area
    overlapArea <- tryCatch(getOverlapArea(refHR, comparedHR), error=function(e) {print(e); return(NA)})

    #' Dahl Overlap Index (2005)
    overlapIndex <- getDahlIndex(comparedArea, refArea, overlapArea)

    #' Percentage of points in polygon
    isPip <- point.in.polygon(comparedLocs[, 1L], comparedLocs[, 2L], 
        refHR@polygons[[1L]]@Polygons[[1L]]@coords[, 1L], 
        refHR@polygons[[1L]]@Polygons[[1L]]@coords[, 2L])
    propInRef <- sum(isPip, na.rm=TRUE) / n * 100

    res <- data.frame(method, percent, n, sizeRatio, overlapIndex, comparedArea,
        refArea, overlapArea, propInRef, stringsAsFactors=FALSE)
    
    return (res)
}


#'#
#'# Compute the overlap area between two surfaces. Returns 0 if no overlap
#'# 
getOverlapArea <- function(s1, s2) {
	intersection  <- gIntersection(s1, s2)
	if (is.null(intersection)){
		return (0)
	} else {
		return(intersection@polygons[[1]]@area)
	}
}


#'#
#'# Dahl overlap Index (2005)
#'# 
getDahlIndex <- function(a1, a2, overlap) {
	return ((overlap / (a1 + a2 - overlap)) * 100)
}


#'#
#'# Create a grid for kernel based on cellsize
#'# 
makeKernelGrid <- function (x, y, cellsize=250L, extent=0.3, extentCell=10L){
    
    minX <- min(x)
    maxX <- max(x)
    minY <- min(y)
    maxY <- max(y)
    
    rangeX <- maxX - minX
    rangeY <- maxY - minY
    
    bufferX <- max(extentCell * cellsize, extent * rangeX)
    bufferY <- max(extentCell * cellsize, extent * rangeY)
    
    Xmin <- minX - bufferX
    Xmax <- maxX + bufferX
    Ymin <- minY - bufferY
    Ymax <- maxY + bufferY
    
    celldimx <- ((Xmax - Xmin) %/% cellsize) + 1
    celldimy <- ((Ymax - Ymin) %/% cellsize) + 1
    
    grid <- GridTopology(c(Xmin, Ymin), c(cellsize, cellsize), 
                         c(celldimx, celldimy))
    spgrid <- SpatialPixels(SpatialPoints(coordinates(grid)))
    
    return (spgrid)
}



prepareHRCompareData <- function(gpsData, argosData) {
    
    res <- vector("list", length(levels(argosData$deployment)))
    for (i in seq_along(levels(argosData$deployment))){
        depl <- levels(argosData$deployment)[i]

        tmparg <- argosData[deployment == i]
        
        if(nrow(tmparg) < 5){
            res[[i]] <- list(ref=NULL, compared=NULL, deployment=depl)
            names(res)[i] <- depl
            next
        }
        
        tmpgps <- gpsData[deployment == i]
        tmpgps <- tmpgps[time < max(tmparg$time) & time > min(tmparg$time)]
        
        spGps <- SpatialPoints(tmpgps[, list(x=longitude, y=latitude)])
        spArgos <- SpatialPoints(tmparg[, list(x=argLong, y=argLat)])
        
        res[[i]] <- list(ref=spGps, compared=spArgos, deployment=depl)
        names(res)[i] <- depl
    }
    return (res)
}


compileHRStats <- function(data, statsCols, statBy) {
    res <- data[, c(lapply(.SD, mean, na.rm=TRUE), lapply(.SD, sd, na.rm=TRUE)), by=statBy,
        .SDcols=statsCols]
    nend <- ncol(res)
    nstart <- nend - 2*(length(statsCols)) + 1
    setnames(res, nstart:nend, c(paste(statsCols, "mean", sep='.'), paste(statsCols, "sd", sep='.')))
    return (res)   
}

